---
title: "R Notebook"
output: html_notebook
---

# Computation of Asymptotic Critical Values
Note: The definition of `mu` here is different from that in the paper. Since we have set variances of each treatment to 1 and sampling rule to 0.5, that means that sigma = 2. The value of mu here should be understood as taking the original mean difference and dividing it by sigma. 

```{r}
#Group sequential test
library(mvtnorm)
library(stats)
library(tidyverse)

conditional_prob_2 = function(gamma, alpha, means, co, mu){
num = pmvnorm(lower = c(-2.797, gamma),
            upper = c(2.797, Inf),
            mean = means,
            sigma = co)[1]
den = pnorm(2.797, mean = mu) - pnorm(-2.797, mean = mu)
return(((num/den) - alpha))
}


gamma = function(mu, alpha){
  #critical_value for first stage
  prob_continuing = pnorm(2.797, mean = mu) - pnorm(-2.797, mean = mu)
  gamma_1 = uniroot(function(x) 1 - pnorm(x, mean = mu) - alpha*(1- prob_continuing),
                  c(2,100),
                  tol = 0.000001)$root
  
  #critical_value for second stage
  means = c(mu,2*mu)
  co = matrix(c(1,1,1,2), nrow=2, byrow=T)
  gamma_2 = uniroot(function(x) conditional_prob_2(x, alpha, means, co, mu),
                  c(1,100),
                  tol = 0.000001)$root
  
  return(c(gamma_1, gamma_2))
}

gamma = Vectorize(gamma)
```


# Plot of Asymptotic Critical Calues
```{r}
#Power envelopes
mu_seq = seq(0, 5, by = 0.5)
df = data.frame(mu_seq = mu_seq)
df = df %>%
  mutate(stage1 = gamma(mu_seq, 0.05)[1,],
         stage2 = gamma(mu_seq, 0.05)[2,]) 
names(df) = c("mu_seq", "First stage", "Second stage")

df = df %>%
  pivot_longer(-mu_seq,
               names_to = "Stage",
               values_to = "Thresholds") 

ggplot(df, aes(x = 2*mu_seq, y = Thresholds, color = Stage)) + #multiply by 2 to get back original mu-value
  geom_point() +
  geom_line() +
  xlab((expression(paste("Drifting null value (", eta[0], ")") ))) +
  ylab((expression(paste("Asymptotic critical values (", gamma, ")" )))) +
  theme(#aspect.ratio = 4/3,
        legend.title=element_blank(),
        text = element_text(size=13),
        axis.text = element_text(size=12),
        axis.title = element_text(size=13),
        legend.position = c(0.85, 0.2)
        )

#setwd("/Users/akarun/Library/CloudStorage/Dropbox/Optimal sequential tests/Figures")
ggsave("Critical_value_thresholds_GST.png")
```

# Monte-Carlo simulation
```{r}
#Simulation with uniform errors

test_simulation = function(n, mu, alpha){
  thresholds = gamma(mu, alpha)
  
  y0 = sqrt(3)*runif(n, -1, 1)
  y1 = (2*mu/sqrt(n)) + sqrt(3)*runif(n, -1, 1)  #multiply by 2 as this is the actual treatment effect due to the scaling of mu described earlier

  sigma0 = sd(y0[1:(n/2)])
  sigma1 = sd(y1[1:(n/2)])
  sigma = sqrt(2*(sigma0^2 + sigma1^2))

  x1 = (sqrt(n)/sigma)*(mean(y1[1:(n/2)]) - mean(y0[1:(n/2)]))
  x2 = (2*sqrt(n)/sigma)*(mean(y1) - mean(y0))

  if (abs(x1) >= 2.797){
    reject = (x1 >= thresholds[1])
  } else{
  reject = (x2 > thresholds[2])
  }

return(reject)
}

finite_n_size = function(n, mu, alpha, M){
  samples = replicate(M, expr = test_simulation(n,mu,alpha))
  return(mean(samples))
}

test_simulation = Vectorize(test_simulation)
finite_n_size = Vectorize(finite_n_size)
```

# Plotting Monte-Carlo for size
```{r}
mu_seq = seq(0, 5, by = 0.5)
M = 25000 

size.df = data.frame(mu_seq = mu_seq)
size.df = size.df %>%
  mutate(size_n50 = finite_n_size(50, mu_seq, 0.05, M),
         size_n100 = finite_n_size(100, mu_seq, 0.05, M),
         size_n200 = finite_n_size(200, mu_seq, 0.05, M),
         size_n500 = finite_n_size(500, mu_seq, 0.05, M))
```

```{r}
names(size.df) = c("mu_seq", "50", "100", "200", "500")

size.df2 = size.df %>% 
  pivot_longer(-mu_seq,
               names_to = "n",
               values_to = "Size") %>%
  mutate(n = factor(n, levels = c("50", "100", "200", "500")))

ggplot(size.df2, aes(x = 2*mu_seq,
                      y = Size,
                      color = n)) +
  geom_point() +
  geom_line() +
  xlab((expression(paste("Drifting null value (", eta[0], ")")) )) +
  ylab("Type I error") +
  geom_hline(yintercept = 0.05, 
             linetype = "dashed",
             size = 0.3,
             color = "blue") +
  theme(#aspect.ratio = 4/3,
        text = element_text(size=13),
        axis.text = element_text(size=12),
        axis.title = element_text(size=13),
        legend.position = c(0.12, 0.78)
        ) 
  #scale_color_brewer(palette = "RdGy")

#setwd("/Users/akarun/Library/CloudStorage/Dropbox/Optimal sequential tests/Figures")
ggsave("Size_GST.png")
```